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Abstract 


In an earlier paper ^ an algebraic characterization 
vas made of the problem of resolving closely spaced 
plane waves incident on a linear array. The character¬ 
ization suggests several data-adaptive processing meth¬ 
ods and encompasses the Wiener, Maximum Likelihood, and 
Pisarenko methods. In this paper, the algebraic ap¬ 
proach is amplified and the results extended to consider 
correlated noise. A recursive algorithm is given for a 
particularly effective processing method. 


An important array processing problem is that of 
determining the directions of propagation of plane waves 
incident on a linear array of uniformly spaced sensors . 2 
Contemporary spectral analysis has led to the develop¬ 
ment of several array processing methods that are able 
to resolve plane waves with nearly identical directions 
of propagation. These methods include the Wiener Pre¬ 
diction Filter method,3 the Maximum Likelihood method,3 
and the Pisarenko method. 1 * This paper amplifies and 
extends an algebraic approach given earlier 1 based upon 
an algebraic characterization of the array processing 
problem. The results encompass the methods mentioned 
above and include the case of correlated noise. A re¬ 
cursive algorithm is presented for implementation of a 
particularly effective processing method. 

Model of the Array Data 

Consider the complex sinsoidal time-space plane 
wave T{t,r) as represented by 

f(t,r) = Ae - U) 

where A is the complex amplitude, t is the continuous 
time variable, r * xi+yj+zk is the continuous space^ 
variable, (v is the (temporal) frequency, and k^ » k x i 
+k J+k k is the wavenumber (spatial frequency). This 
wave travels in the direction of -k with a speed of 

propagation c = yjy . Let us now monitor this wave with 

a linear array placed along the x-axis whereby y ■ z ■ 0 
as shown in Figure 1. The detected signal is 


J[wt+k x] 

f x (t,x) * Ae . (2) 

From this ideal data it is theoretically possible to 
determine the values of the parameters u and k x - Fur¬ 
thermore, if the speed of propagation is a known con¬ 
stant or a known function of frequency, we can then 
calculate the wavenumber's magnitude from 

W • (3) 

Because we do not have com[lete knowledge of the wave¬ 
number it, we cannot unambiguously determine the direc¬ 
tion of propagation. However, we can determine the 
polar angle y associated with the wave as 



This angle defines a cone whose central axis lies along 
the linear array. This information alone is sufficient 
for many applications. For example, the wave may be 
known a priori to be traveling in the xy plane. This is 
the case we shall consider hereafte csr - _; . — 
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Figure 1. Linear array alotl 


We have reduced the problem from four dimensions ir, 
the variables to two dimensions, through the constraint 
of a linear array. We can further reduce the problem to 
one dimension by noting that time and space are inde¬ 
pendent quantities in this model. Thus we can perform 
our analyses for u and k x separately. Ordinary- time 
series processing such as filtering or spectral estima¬ 
tion can be applied to each sensor to first determine 
the presence of signals at a particular temporal fre¬ 
quency u. These outputs, one for each sensor, are then 
spatially processed to determine the directions of 
sources radiating at the frequence u. Hereafter, we 
shall suppress the time domain and consider only the 
3 5 

spatial dimension. 

Figure 2 depicts a linear array of p sensors uni¬ 
formly spaced d units apart. A plane wave is impinging 
upon the array with an incident angle 8. Hoting that 
the incident angle is complementary to the polar angle, 
we have 

k k k 

sin 8 - T-2r = « TT 1 

|k| w/c 

from which it is seen that 

k . ^ - B -f n - 6 (5) 

x A 

where A is the wavelength of the plane wave. Defining 
our origin at sensor zero, the nbh sensor will sample 
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the wave at the point x » nd. Hence, at any particular 
instant in tine the array output ia, from (2) 


yin) 

where t ia a 

in:'* .'lilt. 


■ Ae L * J (6) 

o*n«p 

phase angle dependent on the sampling 


f x (nd) 



Figure 2. Plane wave with incident angle 0. 


"tie set of p instantaneous spatial samples (6) as 
measured by the array is referred to as a "snapshot." 

Ir. this case the snapshot is a sampled complex sinusoid 
whore sampled spatial frequency is given by 

- ■ ? : J ?i n 9 , Clearly, an estimate of the sinusoid's 
frequency directly yields an estimate of the direction 
of rrepagation of the plane wave if the wavelength X is 
Kr. wu. As such, spectral estimation is seen to play a 
imminent role in linear array processing. 

We now generalize our model to include multiple 
p i -it.• waves incident on an array in which the sensor 
in i i ."it ions are contaminated by white measurement noise. 
I f t here are a total of q plane waves and the k** 1 plane 
wave ha. a direction of propagation 0^, it follows by 
.'up rpnr.ition that the snapshot will have the form 


\ " a Jhu. 

y(n) « n(n) + y A^e e , 0£n<_p-l 

k=l 

where the q cinuosoid frequencies are given by 
2nd sin 0. 


(7) 


and n(n) are uncorrelated zero mean random variables 

n 

with variance o‘ that represent the measurement noise. 

We assume that the w, are all different. 

Our objective is to estimate the frequency param¬ 
eters ui using these snapshot measurements. We are 
particularly interested in the ability to resolve, or 
distinguish between two plane waves with very similar 
frequencies (i.e., w i w ^). This estimation capability 
require.-, the utilization of a number of snapshots taken 
sequentially in time. Our data then has the form 


( n) = n ( 


n_ln) 
Pi 


A e 
k 


J *km Jn<1, k 


1 ^m £M, 0 <_n <_ p-1 


k=l 


( 8 ) 


where m ir tin- snapshot index and M i3 the total number 
of snapshots. In this model, we assume that the phase 


angles are uncorrelated random variables uniformly 
distributed on [-*, + *]. Their random nature arises 
from the Independence of the sinusoidal sources and from 
the approximate randomness of tlae-sampllng far below 
the Nyqulst rate. 

It la convenient at this point to represent the 
problem in vector notation. We represent the mth snap¬ 
shot (8) by the p*l column vector 

y * ty (0) y (1) ... y (p-1))' (9) 

"in m in rn 

and define the pure complex sinusoid vector as 

S - [1 e j “ e J? “ ... e J(p ' 1) “]' . (10) 

*Tl) 

Lastly, the noise vector associated with the mth snap¬ 
shot is defined as 

th “ t B (l) ••• h B (p-l)]' • (11) 

With the above notation, we may compactly represent the 
snapshots by the data vector equation 
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km, 


S 
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, 1 <^m <_K 


k-1 


( 12 ) 


The array data (12) is random due to its dependency 
on the random phase angles and the contaminating 
noise n^n). Assuming that these random variables are 
pairwise uncorrelated and statistically invariant with 
respect to the snapshot index m, it follows that each 
data vector v is a vindoved realization of a vide-sense 
stationary raKdom process. The mean value of this proc¬ 
ess is the zero vector. 


UzJ * E^} + 
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Jf 


A.S E(e 
‘"“k 


km. 


- 0 , 

vhile its covariance matrix is specified by 


(13) 
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where 1^ is the p*p identity matrix and * lA^!' 

the power of the kth incident, plane wave. Since th, 
random vector process is vide-sense stationary, the 
covariance matrix R must be positive ser.i-definite, 
Toeplitz, and Hermitian. 

We now describe three contemporary array process ir..- 
methods and then present an algebraic approach to iden¬ 
tifying the frequencies (w^), based upon the structure 
of the data ^ and the associated covariance natrix R. 

Contemporary Processing Methods 


Wiener Filter Method 


The Wiener Filter method is bar.rd on filtering the 
data such that the signal-to-noise ratio at the filter 
output is maximized.8 it is essentially a linear pre¬ 
diction approach that is quite sir.ilar^to the Maxima- 
Entropy method of spectral estimation.' Many adaptive 
array processing algorithms are equivalent to the Wiener 
Filter method, including Alan's orthonom.il lattice 
filter algorithm.^ For the array processing problem an 
optimum weighting vector is obtained by*' 


a = uK'-'h (lb) 

wliere h*[l 0 0 ... 01 ’ and u is nil arbitrary 

scalar. As in the Maximum gntropv Method, the power 
spectrum may be computed by 











p w ( “) 




S + a a + S 


(16) 


Maximum Likelihood Method 


The Maximum Likelihood method is based on filtering 
the data such that power at the frequency of interest!a 
passed and all other frequency components are rejected 
in an optimal manner. In our notation, the power spec¬ 
trum is given by 6*9,10 


P ML W 


s + s _1 s 

—0) “U 


(17) 


Pisarenko Method 


The Pisarenko method has not found as widespread 
use as the Wiener Filter method and the Maximum Likeli¬ 
hood method. Haykin" recently applied part of the 
Pisarenko method-1 to the array processing problem via 
a special autoregressive-moving average (ARMA) model. 

The complete^! sarenko method is based on a theorem of 
Caratheodory that allows decomposition of the exact 
truncated covariance sequence r(n), O^n^q-l, into a 
positively-weighted sum of q complex sinusoids and white 
noise. The method has three steps: 

(i) identifying and removing the noise contri¬ 
bution to the covariance matrix, 

(ii) forming the qxq covariance matrix Rq and 

and analysing the single eigenvector corre- 
pondlng to the unique minimum eigenvalue 

1 . * 0 to determine the sinusoid frequencies, 

min 

and 

(iii) solving a set of q simultaneous linear 
equations for the sinusoid powers. 

Algebraic Processing Approach 

We now formulate a generalized minimization problem 
which suggests several particular methods. Under dif¬ 
ferent constraints, the solution of this problem encom¬ 
passes each of the methods in the previous section. 

Let us consider a general nontrivial p*l coeffi¬ 
cient vector a that is orthogonal to the noise-free 
component of each of the data vectors From (12), 

this orthogonality is defined by the inner product 
relationship 

0 * 

■J]V <■• 8, > U8) 

k«l k 

Since the (u. } are all different and the (♦, ) are ran- 
K Kin 

dom in nature, a little thought will convince oneself 

that a must be orthogonal to each of the q sinusoid 
vectors S , 1< k< q. We next define the general z- 
^k 

transform A(z) of the coefficient vector a by 

A(z) » <a, z»> (19) 

where z * [l z* 1 z _? ... z 1_p ] . It is then readily 

shown that the orthogonality of a to each S , 1 < k < q, 

”**k 

implies that A(z) must have q finite zeros located on 

J “k 

the unit circle at the points z^ = e , l^k^q. With 
this in mind, the required frequencies can be deter¬ 
mined by examination of the zeros of A(z). 


In the Idealized noise-free case, the snapshot vec¬ 
tors are members of the q-dimensional subspace span¬ 
ned by the q linearly independent sinusoid vectors S , 

”®k 

l£k£q. Thus, forq<pthere always exists a p*l vector 
a that is orthogonal to the noise-free components of the 
snapshot vectors y^ . If q>_p, there generally does not 
exist such a vector a. Furthermore, when noise is 
present, even if q<p, there generally doe: not exist a 
vector a that is orthogonal to each of the noise- 
contaminated snapshot vectors jfa. nonetheless, it is 
intuitively desirable to select a coefficient vector 
which is nearly orthogonal to each of the snapshot vec¬ 
tors in some well-defined manner, and to determine the 
plane wave sinusoidal frequencies by examination of the 
zeros of the z-transform of this coefficient vector. A 
convenient method to evaluate these zero locations is to 
search for nulls in the magnitude of the Fourier trans¬ 
form of the coefficient vector. Since there can be more 
zeros than plane waves (i.e., p-l>q), we can estimate 
the (P^) in order to separate "signal zeros" from "noise 
zeros". 

To obtain a mathematical measure of closeness to 
orthogonality, it is beneficial to define an orthogo¬ 
nality error vector e(a) whose mth element is the inner 
product of a and y^ . denoted by < a,y r > . We define the 
optimum a to be the vector a° that minimizes some posi¬ 
tive definite functional f of £(a). Hence we write 

e(a) = [e(l) e(2) ... e(M)] 1 

where 

e(m) = < a,^ > 

and 

f(e(a°)) * min f(e(a)) 
atA 

where A is a constraint set for the solution vector a°. 

A constraint is generally necessary for the minimization 
to be well-defined, that is, for a° to be unique and 
nontrivial. 

We must choose an inner product for (20) and ar. 
error functional f for (21). Let us choose in particu¬ 
lar the standard vector inner product < a,^ > * a'y in 
which case we have 


( 20 ) 

( 21 ) 



( 22 ) 



A convenient positive definite functional for an error 
vector is the mean square error criterion 

f(e) • E{ || e || P ) (23) 

where || e || » ^|e(l)| P ♦ ... ♦ |e(M)| P , the Euclidean 

norm of e.l** It will be computationally expedient to 
normalize this criterion by the length M of the vector 
e. From (22) and (23) we have 



■ a*Rn (2l) 
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where R is the covariance matrix defined in (lit). Apply¬ 
ing (21), the quadratic form (2it) must now be minimized 
according to some constraint that causes t»° to be unique 
and nontrivial. Next we consider two possible con¬ 
straints, a linear constraint and a quadratic con¬ 
straint. 

Linear Constraint 

The first constraint is that a 0 lies on a hyper- 
planc specified by 

A • (at C P :a + h+h t a » 2) (25) 

where h is a nontrivial pxl vector that characterizes 
the orientation of the hyperplane. The solution to 
(21) with this constraint is given by 


1° *r-rVl R ' 1 ii 

h R h 


and the criterion's minimum value is 

f(a°) « * . (27) 

h R h 

Quadratic Constraint 

The second constraint is that a° lies on the quad¬ 
ratic surface specified by 

A * (a< C P :a + Wa - 1) (28) 

where W is a positive definite, Hermitian matrix which 
characterizes the quadratic surface. The solution to 
(21) with this constraint is given by 


/x + . Wx . 
i -min -min 


and the criterion's minimum value is 


Extension to Correlated Noise 

Through selection of the constraint set, the alge¬ 
braic approach extends quite readily to the case in 
which the contaminating noise is correlated. Such is 
the case when the noise is due not only to sensor meas¬ 
urement noise but also to a directional background 
noise field in the array environment. Note that any 
undesired signal (jasning interference, for example) 
may be considered as correlated noise. 

Generalizing (lL), we have that the data covariance 
matrix for the case of correlated noise is given by 


R * o B ♦ 


s + 

-v-k 


where the noise covariance matrix B is defined by 

o 2 B • E(n r/) . (32) 

TB m tK 

We assume that the shape of the noise spectrum is known, 
which implies knowledge of B. 

Returning to the quadratic constraint (29,30), we 

note (X . , x . ) is the solution to 
min* -min 

(R-X . W)x , * 0 . (33) 

min -min — 

We have seen that for the choice W = I we have the 

P 

Pisarenko solution. In this case it is well known that 

2 

X , =o and that we essentially have a white-noise 
min -« 

power cancellation algorithm. Hence it is a simple 
step to choose 


. . and achieve a colored-noise correlation cancellation 

algorithm. This step can be Justified further by- 


algorithm. This s 
rewriting (31) as 


(R-o 2 b) 


S S + 
^V“k 


where (X ,x , ) is the minimum-eigenvalue and eigen- 
min -min 

vector pair of W^F. 

The above solutions encompass the three contempo¬ 
rary methods described earlier. Forh * [l 0 ... 0]', 
Z;uation (26) is precisely the Wiener Filter solution 


h + R -1 h 


Just as in linear prediction, this 


constraint implies that the first element of a 0 is 
fixed at 1 and the other elements are unconstrained. 

For h = F , (27) is precisely the Maximum Likelihood 
solution. This constraint requires A°(z) to have unity 
gain at z * eJ“, while the minimization strategy opti¬ 
mally reduces the gain at other frequencies. For 
W = Ip, the quadratic surface is the hypersphere of 
radius one, and ( 29 ) is a generalized version of the 
Pisarenko method. There are several differences. 

First, no special ARMA model is invoked, as is done by 
Haykin. 1 ' Second, neither noise power removal nor ma¬ 
trix order reduction are required, as they are In the 
Pisarenko method. Third, this method is based upon a 
minimization strategy and so Justifies estimates, gen¬ 
erally even non-Toeplitz, of the covariance matrix R. 

In the special case of a Toeplitz estimate matrix, a 
power identification technique similar to the Pisarenko 
method ear be employed, as is shown later. Fourth, the 
general constraint matrix W allows greater flexibility 
in the solution. The quadratic constraint solution 
generalizes the Pisarenko method and extends it to the 
multipie-snapshot array processing problem. 


and remembering that we seek a vector that is orthogonal 

to the sinusoid vectors (S ). Also, it is apparent 
'“k 

+ * 

from the constraint a Ba = 1 that we are again specify¬ 
ing a set of vectors with constant norm, but now the 
norm is determined by the noise covariance matrix B. 

The linear constraint (26,27) may also be extended 
for correlated noise. We shall only consider the 
Wiener solution, since it has been shown to achieve _ 
greater resolution than the Maximum Likel ihood Method. ’ 
We will show later that to extend the Wiener linear pre¬ 
diction solution, a reasonable constraint h given B is 


Note that we no longer have a linear predictor. 

Example 

To illustrate the linear and quadratic constraint 
solutions, let us consider the case of a single plane 
wave of power F. and spatial frequency u incident on 
an array of two sensors in a correlated noise field. 
For this case the noise covariance matrix is given by 







2 2 
o B = a 


[::] 


where |b|<l, and the data vector covariance matrix by 


2 

Pe 1 «• o b* 


'^e 1 ♦ oT) 


First we consider the linear and quadratic con¬ 
straint solutions without accounting for the correlated 
noise (i.e., the Wiener and Pisarenito methods). For 
the linear constraint we choose = [l 0 ... 0]' 

and find from ( 26 ) that 


P e 1 + o^b 


Pefir.irg the signal-to-noise ratio by SHR 
see that A°(z) has a zero located at 


LI SNR + b e 1 

L SKR + 1 J 


>1 2 , 

IP^ 1 ♦ o Z b| 


Thus A°(z) has a zero located at 


OT e + b 

Z * J 

Ju. 


We aee in this case that the zero lies directly on the 
unit circle, regardless of the signal-to-noise ratio. 
In fact, when the noise is white, the zero perfectly 
indicates the plane wave spatial frequency u^. How¬ 
ever, when the noise is correlated, there is a fre¬ 
quency bias present that again becomes greater as the 
signal-to-noise ratio decreases. 

Let us consider now the linear and quadratic con¬ 
straint solutions which account for the correlated 
noise. For the linear constraint we choose 
h » B(l 0 ... 0)' and find that 


Ve A -b) 

P^l-be 1 ) + o(l-|b| 


where ii is a scalar function of P. , u , o , and b. 
Thus A°(z) has a zero located at 1 1 


.1 SNPtl-be *) 

-Jw. p 

SHR(l-b e A ) + l-|b| 


As in (1.0), we see that pure white noise will still 
cause the zero to migrate away from the unit circle, 
and that correlated noise will introduce frequency bias. 
However, as the noise becomes "more correlated" (i.e., 
|b| ■*1), the zero moves closer to the unit circle and 
asymptotically indicates the exact plane wave spatial 
frequency regardless of the signal-to-noise ratio. 
Hote that the effect of an interfering harmonic source 
J“>p 

(i.e., b • e ) is completely removed. 

For the quadratic constraint we choos V = B arid 
find that 


As has been noted elsewhere,*5 this linear- 
prediction solution suffers from zero migration away 
from the unit circle even when the noise is white (i.e., 
b = 0). This migration degrades resolution, since 
applying the Fourier transform to evaluate zero loca¬ 
tions may indicate only a single null when in fact 
there are two zeros close together somewhat off the 
unit circle. Furthermore, we see there is a frequency 
bias introduced by the correlated noise. This bias 
becomes greater as the signal-to-noise ratio decreases. 

For the quadratic constraint we choose W = I and 
find 'Yom (29) that ** 


J“l 2 
P,e + o b 


where v Is a Bcalar function of and b. Thus A°(z) 
has a zero located at 


We see that the zero indicates the exact plane wave 
spatial frequency « , regardless of the signal-to-noise 
ratio or the particular value of b. For this reason, 
we expect the quadratic-constraint solution to obtain 
high resolution. 

To summarize the development to this point, the 
algebraic approach is based on approximating an orthog¬ 
onality condition between a solution vector and each 
of the data vectors. This approach encompasses three 
contemporary array processing methods and readily 
extends to the case of correlated noise. 


Implementation of the 


dratic-Constralnt Solution 


In the previous section we saw that the quadratic- 
constraint solution (29,30) is a promising array proc¬ 
essing method in terms of its perfect resolution given 
exact covariance values. However, it requires an 
eigenvalue-eigenvector computation that seems to be 
quite burdensome. Fortunately, a simple recursive algo¬ 
rithm can be derived using, the nature of the array 
processing problem. , 

First we recall the standard "inverse iteration" 11 ' 
method for finding the minimum eigenvalue and eigen¬ 
vector pair of a complex matrix D. Consider the se¬ 
quence of vectors (x^) defined by 

D *w ‘ \-i (1,T) 

where Xq is a nonzero and arbitrary. Ac K increases, 
we have 


5 









It can be shown that 6X Is approximately given by 


^k Tain 


and 


Vk-1 

t 


*min 


This method is appropriate to the array processing prob¬ 
lem, in which the data arrives sequentially. Assume 
that from M snapshots we have estimated the covariance 

matrix by R^ and obtained the desired pair (^, ^ 1, 

When the next snapshot is available, we form R M+ ^ and 
compute from (1*7) using D * R M+1 and x^ » x^. 

Since the inverse iteration method generally has fast 
convergence, a single Iteration of (1*7) for may t> e 

sufficient as long as is only slowly time-varying. 

To accelerate convergence, we can apply the "inverse 
iteration of Wielandt" 1 ^ wherein an approximation of 
*min is subtracted from the main diagonal of D before 

iterating. Given R„ ., we use X to approximate X u . 

M+l m 1 

The iteration is given by 


(R M+1 ■ X M I)5 M+1 * ’ 

3m* 1 ■* *nin * and 


^I+A+l 


= A -*■ A . . 

M+l min 


(1*8) 


For a Toeplitz and Hermitian covariance matrix estimate, 

2 

each iteration can be performed with 0 (p ) multiply- 
adis using Zohar's algorithm. 17 An alternate algorithm 
by Oueguen 1 ® can be used to avoid numerical difficulties 


that 


may be associated with (R M+1 - * M I). 

For the case of correlated noise, we have 


D 


r-i: 


' M+l ’ 


iteration to 
iteration is 


We suggest generalizing the accelerated 

avoid calculation of P Our general 
given by 


(R M+1 ~ V 'b+l ^ = % 


M ’ 


4l + l " 4un > 

M t 


M*1 


and 


min 


(1*9) 


. x (4D)x 

«> ■ --- (r ? ) 

xWx 

Our approximation to the new X is then given by the 

sum of 4X and the previous X , . With appropriate 

Bln 

definitions, this approximation replaces X in ( 1 * 9 ) 
when we expect the new covariance matrix estimate to 
differ considerably from the previous estimate. 

Relationship to Linear-Constraint Solution 

The iterative implementation of the quadratic- 
constraint solution gives insight into the linear- 
constraint solution. Namely, the first iteration of 
(1*7) with D»R and x^ * h yields the linear-constraint 
solution of ( 26 ) within a constant of proportionality. 
Repeated calculation of the linear-constraint solution, 
with h at each step equal to a° of the previous step, 
is in fact an iterative implementation of the quadratic- 
constraint solution. It is apparent that at each step, 
the constraining hyperplane is realigned according to 
the estimated solution. With these Insights, it is 
reasonable to choose 



(53) 


as the linear constraint for correlated noise, since it 
yields the first step of the iterative quadratic- 
constraint solution (without acceleration) for corre¬ 
lated noise given in (1*9). This Justifies the choice 
made earlier in ( 36 ). 

In this section we have presented a recursive 
algorithm (1*9) for implementing the quadratic-constraint 
solution. The algorithm makes use of the sequential 
nature of the snapshot data to efficiently employ in¬ 
verse Iteration. The algorithm includes the case of 
correlated noise. A modification to the algorithm (52) 
was presented for the case where successive covariance 
matrix estimates differ considerably. 

Covariance Matrix Estimate 


2 

Each iteration still only requires 0(p ) multiply-adds, 

and B is never calculated. 

In some applications, it may not be desirable to 
calculate the solution vector after every snapshot. 

For instance, forming a new covariance matrix estimate 
and calculating a new solution only every L snapshots 
reduces the average computation rate by a factor of L. 
Unfortunately, if the new covariance matrix estimate 
differs considerably from the previous, the previous 
eigenvalue may be a poor approximation to * m j n * and 

convergence will be slowed. To obtain a better eigen¬ 
value approximation, we apply perturbation techniques.^ 
Suppose that 0 and W are Hermitian matrices and that we 
have solved the eigenvalue-eigenvector problem 

Dx * XWx . (50) 

Applying a Hermitian perturbation 40 to D we have the 
new problem 

(D+4D)(x+4x) * (X + 4X )W(x ♦ 4x) . (51) 


To employ the proposed processing methods, an 
estimate of the covariance matrix is required. From 
this estimate, a solution vector i ; obtained and the 
zeros of the vector's z-transform examined to determine 
the plane wave spatial frequencies. Given a p*l solu¬ 
tion vector and q plane waves, q<p, there will be q 
"signal" zeros and p-q-1 "noise" zeros. These zeros 
must be separated from one another. It is well known 
that in the linear prediction solution, dominant fre¬ 
quency components will generate zeros closer to the unit 
circle than less powerful components; thus, a simple way 
to evaluate signal zero locations is to search for nulls 
in the solution vector's Fourier transform. For the 
quadratic solution in white noise, it can be shown^^ 
that all of the zeros will be on the unit circle when 
the covariance matrix estimate is both Hermitian and 
Toeplitz. Thus the estimated frequencies can be di¬ 
rectly employed in a power determination technique' 1 * ' 
and the zeros separated on a basis of signal power as 
before. 


6 







3 

A standard covariance matrix estimate is 
M 

vsEyl- (53) 

m=l 

This estimate is unbiased. Hermition, but in general 
not Toeplitz. Furthermore, only one lag product from 
each data vector is used in formulating each element of 
R M - An alternate estimate is the matrix R^ whose ele¬ 
ments are given by 

R M (i,j) «= c(i-j), l<i,j<p (5M 

where 

M p-n-1 

c(n) * m]C^ 2 y a (t+n)y I (tK °i"- p - 1 

m=l 1=0 

c(n) = c # (-n) , -p + l±n<0. 

Tills estimate is unbiased, Hermitian, and Toeplitz. 
Also, p-n lag products from each data vector are used 
in formulating each element c(n). Thus the estimate of 
(5 1 *) has lower variance than that of (53). 

In the following simulations, the standard non- 
Toeplitz estimate 1^(53) will be used in the linear- 
constraint solution in order to compare with previous 
simulations.3 However, for the quadratic-constraint 
solution the Toeplitz structure is important, hence the 

Toeplitz estimate R (51*) will be used. 

M 

Simulation Results 

To compare the performance of these two processing 
methods, the data vectors (12) were generated by com¬ 
puter simulation. The simulation model corresponded to 
that chosen by Gabriel^ in his comparative paper. 

Namely, the case of two sources incident on an array 
with white noise was considered. The parameter selec¬ 
tions were q= 2, p = 8, a 2 = 1, = A ? = 31-62 (30 dB SNR) 

and 3.162 (10 dB SNR), 6^18°, 6 2 = 22°, M=50 (many 
snapshots) and 10 (few snapshots), and d = X/2. 

With vhite noise, the linear solution and the 
quadratic solution correspond to the Wiener and 
Pisarenko methods. The simulation results are shown in 
Figure 3. In this figure, the linear solution has been 
evaluated via its Fourier transform and the quadratic 
solution via the power determination technique. Over- 
layed solutions for ten different realizations of the 
random data are shown to give a sense of each method's 
consistency. 

These results show that both methods work well at 
high SNR with many snapshots. However, the linear 
solution performs poorly at low SNR with few snapshots, 
while the quadratic solution continues to give good 
resolution and good suppression of spurious effects. 

In general, the quadratic solution has shovn better 
performance than the linear solution over a vide range 
of conditions.^ Further simulations are underway to 
compare the performance of the methods in correlated 
noise and to evaluate the recursive algorithm presented 
above. The results will be given at the conference. 

Conclusions 

We have detailed an algebraic approach to array 
processing based upon approximation of an orthogonality 
condition. The approach encompasses several contempo¬ 
rary, high resolution methods. Previous results were 
extended to the case of correlated noise, and a recur¬ 
sive algorithm presented for the quadratic-constraint 


solution. The quadratic-constraint solution appears to 
he particularly effective and suggests further investi¬ 
gation of eigen-analysis array processing methods and 
their implementation. 
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Figure 3. Two-source simulation with sources at 18 and 22 degrees. 

(A) 30 dB SNR, 50 snapshots | 

(B) 10 dB SNR, 10 snapshots / Linear solution 

(C) Single trial from (B) ) 

(0) 30 dB SNR, 50 snapshots } 

(E) 10 dB SNR, 10 snapshots 

(F) Single trial from (E) 


Quadratic solution 
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